Role of quantum nuclei and local fields in the x-ray absorption spectra of water and ice 
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We calculate the x-ray absorption spectra of liquid water at ambient conditions and of hexagonal 
ice close to melting, using a static GW approach that includes approximately local field effects. 
Quantum dynamics of the nuclei is taken into account by averaging the absorption cross section 
over molecular configurations generated by path integral simulations. We find that inclusion of 
quantum disorder is essential to bring the calculated spectra in close agreement with experiment. In 
particular, the intensity of the pre-edge feature, a spectral signature of broken and distorted hydrogen 
bonds, is accurately reproduced, in water and ice, only when quantum nuclei are considered. The 
effect of the local fields is less important but non negligible, particularly in ice. 
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In the last decades high resolution core level spec- 
troscopy, such as x-ray absorption spectroscopy (XAS) or 
x-ray Raman scattering (XRS), has emerged as a power- 
ful tool for the investigation of matter at the molecular 
scale Jj]. Important studies of liquid water and other 
disordered hydrogen-bond (H-bond) environments have 
been performed with these techniques [2|-|^. Core spec- 
troscopy is an element-specific local probe, complemen- 
tary to diffraction for structural information, but requires 
a good theory of the excitation process as aprerequisite 
to a consistent experimental interpretation @1 ■ Given 
that disordered environments need large cells and aver- 
ages over different realizations, computationally efficient 
schemes based on density functional theory (DFT) have 
been the technique of choice in this context Q- DFT, 
however, is not an excitation theory and yields spectra 
that are at best in semi-quantitative agreement with ex- 
periment. 

Dramatic improvement in the calculated spectra is ob- 
tained with a proper excitation theory, as shown in a 
recent paper [8|, which focused on water and adopted 
Hedin's GW approximation 9] for the self-energy of the 
excited electron in presence of a static core hole. Com- 
putationally, this approach is significantly more demand- 
ing than DFT and Ref. H compromised on the structural 
model and the screening approximation. The calculated 
spectrum deviated from experiment mostly in the pre- 
edge intensity which was underestimated by ~50%. It 
was unclear if this discrepancy was due to limitations 
of the structural model, the screening model, or both. 
Ref. ^ used liquid structures generated by ah initio molec- 
ular dynamics (AIMD) in a periodic cell with 32 H2O 
molecules. With standard DFT approximations and clas- 
sical nuclear dynamics this scheme yields enhanced water 
structure. Thus the temperature of the simulation was 
set to T ~ 360i^ to improve agreement with experiment 
at ambient conditions (T — iOQK). In addition, a simple 
uniform screening model was adopted, neglecting the mi- 



croscopic inhomogeneity of the medium. The same paper 
reported also the spectrum of a proton ordered cubic ice 
structure that minimized the DFT energy at T = QK . 
The deviation from experiment was significantly larger 
than in water, but one should recall that in real cubic 
(Ic) or hexagonal (Ih) ices, the distribution of the pro- 
tons in the H-bonds is disordered. Moreover, the effect 
of zero-point motion was ignored. This is large, reflect- 
ing the quantum character of the nuclei, particularly the 
protons, signaled by neutron scattering experiments 

In this paper we address the above issues. We take 
the quantum character of the nuclei into account us- 
ing molecular structures generated by path integral ah 
initio molecular dynamics (PI-AIMD) [lafor water at 
T = 300K and for ice at T = 269K^. In addi- 
tion, we improve the treatment of screening by includ- 
ing, albeit approximately, the local fields due to the in- 
homogeneity of the medium within the Hybertsen-Louie 
(HL) ansatz IJ]. We adopt larger simulation cells than 
in Ref. H and reduce the statistical error in the aver- 
ages. Quantum effects bring the pre-edge intensity in 
close agreement with experiment in both water and ice. 
In addition, they reduce the post-edge intensity, for a 
better agreement with experiment, especially in water. 
Local field effects are modest in water, but cannot be 
neglected in ice that has a more open and less homoge- 
neous microscopic structure. Here they improve the spec- 
tra in the main region that includes near- and post-edge, 
but some residual discrepancies remain, suggesting that a 
more accurate screening model could further improve the 
experimental comparison. Even at this stage, however, 
our approach, combining accurate ah initio simulations 
and spectral calculations, is a powerful interpretive tool 
for x-ray spectroscopy. 

We base our calculations on the molecular structures 
of Ref. ITi, which were obtained with PI-AIMD simula- 
tions using the BLYP functional [l5|, l20| for the elec- 
tronic energy within DFT. The corresponding classical 
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AIMD structures with the same functional approxima- 
tion and thermodynamic conditions were also reported in 
the same paper. Quantum nuclei soften the structure of 
water predicted by DFT, thereby greatly improving the 
experimental agreement of the calculated pair correlation 
functions. Residual overstructuring can be attributed 
to limitations of the functional approximation. Quan- 
tum dynamics increases the fraction of broken bonds and 
broadens the distribution functions. Thus, empirical ad- 
justments like increasing the temperature of the simula- 
tion as in Ref. 18| can be avoided. Ref. |l3| used open PI 
trajectories [3] to calculate the momentum distribution 
of the protons, finding excellent agreement with deep in- 
elastic neutron scattering experiments. Relative to its 
classical counterpart, the quantum momentum distribu- 
tion is strongly displaced towards higher momenta and 
deviates significantly from a spherical Gaussian distribu- 
tion [17[ , indicating that a higher effective temperature in 



classical simulations can only be a poor surrogate of full 
quantum simulations. The latter are based on Feynman 
paths in imaginary time, which provide an exact repre- 
sentation of equilibrium quantum statistical mechanics. 
The x-ray absorption cross section is a dynamic prop- 
erty requiring, in principle, path components in real time 
which would carry a complex phase not amenable to com- 
puter simulation, but we can rely on the wide separation 
of time scales between nuclear and electron dynamics, in 
the spirit of the Franck-Condon principle. In the present 
case it amounts to compute the spectra for the frozen 
nuclear configurations sampled with the adopted discrete 
representation of the Feynman paths in imaginary time. 

Molecular configurations for water at T = 300K and 
ice Ih at T = 269^ were generated in Ref. [13 by adopting 
the "bulk" approximation in which the Feynman path of 
one proton per molecule is left open to improve the statis- 
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was shown 



tics of the momentum distribution 
that this approximation has negligible effect on the mo- 
mentum distribution and very minor effect on the spatial 
correlations between the protons [16[. The simulations 



used a periodic cell with 64 molecules for water and 96 
molecules for ice Ih in a proton disordered configuration 
with cell volumes fixed to experiment. The Feynman 
paths were discretized using 32 replicas in both cases. In 
our XAS calculation we use the static Coulomb hole and 
screened exchange (COHSEX) approach of Ref. and 
average the excitation cross section over all the molecules 



in the cell. We include quantum effects by averaging over 
path-integral replicas. Each particle {H or O nucleus) 
path runs in imaginary time from r = to r ^/ksT. To 
minimize the effect of the open paths we only include in 
the average 26 middle path replicas in water and 16 in ice 
as this was sufficient for convergence. Even limiting the 
average to a fraction of the replicas, the number of cal- 
culations is huge and we consider only a single molecular 
dynamics snapshot to reduce the computational burden.. 
We checked for classical nuclei that this has minor effect 
on the spectra. Aiming at further reducing the compu- 
tational cost, we have investigated whether calculations 
at the centroids of the paths could capture most quan- 
tum effects. Unfortunately, this was not the case. In 
water, centroid calculations led to a pre-edge intensity 
approximately half-way between classical and quantum 
data, consistent with the observation that the fraction 
of broken bonds is approximately the same in the cen- 
troid and path structures. In ice, that has no broken 
bonds, the centroid spectrum was indistinguishable from 
the classical one. 

We include local fields in the screening by means of the 
HL local density approximation according to which the 
screened interaction reads 
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W{v, r') = i {W[r - r'; p{r')] + W[r' - r; p(r)]) (1) 

Eq.(l) is based on the observation that the screening 
strength generally follows the local charge density. Here 

I^[r'-r;p(r)] = ^ J e-i[q; p(r)]t;(q)e''i-("--') rfq 

(2) 

and V is the bare Coulomb interaction. For the dielectric 
function, we adopt Bechstedt model [l8| . 



6[q,p(r)]^l + [(6o-l)-i + a(^)2 



(3) 



where q^ f and LOp denote Thomas-Fermi wavevector and 
plasmon frequency, respectively. These quantities de- 
pend on the local density but in Ref. Q dependence on 
the average density was assumed. The q = Q value eq is 
taken from experiment and a is fixed by requiring Bech- 
stedt model (Eq.(3)) to match the dependence of Penn 
model [l9( . Eq. ^ can be analytically evaluated yielding 



W[v' - r;p(r)] = "^^^^ ^ ^ (4) 

eo a(a;i - a;2)|r' - r| y xi X2 j 

I 



where xi^2 — {—b ± vP — 4ac)/2a and a — l/Aujp, b — a/q^p and c = £0/(^0 — 1)- There are two 
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tributions to the screened interaction W in Eq. (|4]), the 
bare interaction v divided by the macroscopic dielectric 
constant, and a local-density dependent screened inter- 
action, which is fully nonlocal but short-ranged in space. 
To evaluate the effect of W on the self-energy we use the 
convolution theorem for the part depending on the bare 
interaction, which is conveniently calculated in Fourier 
space, and we compute the integral in real space for the 
short-ranged nonlocal contribution. 

Fig.[l]shows the calculated XAS intensity for liquid wa- 
ter and ice Ih at different levels of approximation. The 
spectra are aligned at the onset and the absorption sum 
rule is imposed to fix the normalization so that the spec- 
tra have equal area in the energy range of the figure. 
Hartree-Fock results are also reported: they correspond 
to calculations with the dielectric function set equal to 1, 
i.e. ignoring screening. 

Replacing classical with quantum nuclei affects the 
spectra, because zero-point motion significantly enhances 
the positional disorder relative to classical thermal mo- 
tion. Qualitatively similar conclusions were reached in a 
recent study of the XAS spectrum of ice using PI config- 
urations generated with an empirical intermolecular po- 
tential [20|. In addition, quantum fluctuations increase 
the fraction of broken H-bonds in the liquid by ~4% [3] ■ 
As a consequence pre-edge is increased, post-edge is re- 
duced and there is an overall slight broadening of the 
spectrum. The effects on pre- and post-edge are more 
pronounced in the liquid, where the post-edge is partic- 
ularly interesting as it correlates with the less overstruc- 
tured H-bond network of the quantum simulations. 

Local fields have comparatively smaller effects than 
quantum nuclei. Here it is useful to start from Hartree- 
Fock. In this limit screening is absent, the quasi-particle 
excitation energies are overestimated resulting in too 
wide spectra with the post-edge more prominent than the 
near-edge in both water and ice. COHSEX calculations 
improve the agreement with experiment, by reducing the 
overall width and by lowering post-edge below near-edge 
in the liquid. The difference between the spectra with a 
uniform screening model and those obtained with the HL 
local density approximation is a measure of the inhomo- 
geneity of matter at the molecular scale, which renders 
screening less effective than in the corresponding uniform 
medium. Local fields are underestimated in the HL ap- 
proach and one should expect that screening would 
be further reduced in a more accurate theory. The effect 
should be small in water where the HL induced change is 
minor, but in ice the HL induced change is non negligible 
suggesting that a better screening model should move the 
spectrum further in the Hartree-Fock direction.. 

In Fig. [2] our results with quantum nuclei and COH- 
SEX with local fields, are compared to XAS and XRS 
experiments. The spectra are aligned at the onset and 
normalized by the same area. The difference spectrum 
in the bottom panel of the figure is obtained by sub- 
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FIG. 1: (Color online) Calculated XAS spectra for liq- 
uid water at T = 300K (top) and ice Ih at T = 269K 
(bottom) according to four different schemes: (a) COHSEX 
with classical nuclei and uniform screening (black) ; (b) COH- 
SEX with quantum nuclei and uniform screening (blue); (c) 
COHSEX with quantum nuclei and local field (LF) screen- 
ing (see text) (red); (d) Hartree-Fock with quantum nuclei 
(green) . The spectra are area normalized and Gaussian broad- 
ened(fwhm=0.4 eV). In order of increasing energy, three fea- 
tures, the pre-edge, the near-edge, and the post-edge, charac- 
terize the spectra. 



tracting ice from water intensity dividing the result by 
the integrated spectral area. The accord among experi- 
ments is excellent: the only noticeable difference between 
XAS and XRS in the post-edge of ice may reflect differ- 
ences in the experiments and the samples. Overall the 
agreement between theory and experiment is quantita- 
tively good in both water and ice. In the latter case, 
the dramatic improvement relative to Ref. Q is a strong 
indication of the importance of a realistic model of dis- 
order not only in liquid but also in crystalline water. 
The agreement between theory and experiments in the 
pre-edge is particularly satisfying given that most previ- 
ous calculations failed to reproduce the correct relative 
intensity of this feature, associated to a bound exciton 
mostly localized on the excited molecule [l^ . Based 
on the quantum-classical trends in Fig. [T] we should ex- 
pect that isotope effects should be detectable: making 
the protons more classical should weaken the pre-edge 
and strengthen the post-edge. Remarkably, this has been 
observed in high resolution XAS measurements on liquid 
H2O and D2O samples [i^. The main remaining differ- 
ence between theory and experiment in the water spec- 
trum is a small overestimate of the post-edge intensity. 
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FIG. 2: (Color online) Comparison of calculated XAS spec- 
tra with experiments for water (top), ice Ih (middle) and the 
difference between water and ice Ih (bottom). Two sets of ex- 
perimental data are included for each system [1,0, [13]. Water 
is at room temperature in theory and experiments, ice is at 
T = 2mK in the present calculation, at T = imK in [U 
and at T = 40-?^ in The difference spectrum is obtained 
from [I(water)-I(Ih)]/S where I denotes spectral intensity and 
S is the integrated area of the spectra in the energy range of 
the figure. 



This may reflect residual overstructuring present in the 
H-bond network of the PI-AIMD simulations, a conclu- 
sion also borne by the analysis of the radial distribution 
functions 
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We think that this inaccuracy is due to 
limitations of the BLYP functional approximation. Im- 
proved approximations including hybrid exchange that 
reduces the self-interaction error 2J| and including dis- 
persion (van der Waals) interactions would further soften 
the network structure [2^ weakening the post-edge fea- 
ture. The difference between ice and water in the rela- 
tive strength of near- and post-edge features has its ori- 
gin in the partial collapse of the H-bond network that 
brings a non-bonded molecular fraction within range of 
the first coordination shell in the liquid [8]. The corre- 



sponding "interstitial" molecules do not form H-bonds 
with the central molecule but are H-bonded to molecules 
outside the coordination shell. This interpretation is con- 
sistent with diffraction studies (26j and we confirm it 
here. Quantum nuclei are crucial in water because they 
enhance the network collapse compared to classical sim- 
ulations. This effect is absent in ice Ih, which has an in- 
tact tetrahedral network with coordination 4 27[ . In ice 
the residual differences between theory and experiments 
in the main edge may reflect to some extent the differ- 
ent temperatures of the simulation and the experiments. 
However, thermal effects should be small as the vibra- 
tional spectrum is largely ground-state dominated [it}. 
We think that the differences should mostly reflect the 
adopted screening approximation, which underestimates 
local fields. A more accurate screening model should en- 
hance post-edge at the expense of near-edge, improving 
the agreement with experiment. A fully first-principles 
treatment of static screening would be possible but com- 
putationally very demanding. Finally, we have neglected 
the energy dependence of the self energy, which should 
affect the calculated spectra, particularly at high en- 
ergy. Such effects could be included either via a dynamic 
COHSEX calculation or, simply, with phenomenological 
schemes such as rescaling the energy axis or using an 
energy- dependent broadening function 2^, 2^ . The good 
agreement between experiment and theory found in the 
present study suggests that these effects should be small 
in water and ice. 

In conclusion, we have improved an approach to calcu- 
late the XAS spectra of disordered environments [8]], by 
including local field effects beyond the uniform screening 
approximation. The most important result of the present 
investigation has been, however, to elucidate the role of 
nuclear quantum dynamics on the spectra of water and 
ice. On one hand zero point motion broadens the spec- 
tra via Franck- Condon factors that we calculate by av- 
eraging over the distribution of the Feynman paths. On 
the other hand quantum nuclei modify the liquid struc- 
ture by facilitating H-bond breaking. Our analysis shows 
that many competing factors are important to reproduce 
the spectral differences observed between water and ice. 
In view of its ability to probe local environments, x-ray 
spectroscopy is a very important tool for investigating 
not only bulk neat water, as done here, but also water 
solutions, water at interfaces, and, in general, soft con- 
densed matter environments like bio-molecular solutions. 
The theoretical techniques developed here should con- 
tribute to these studies by giving insight on the origin of 
the spectral features and by allowing us to connect these 
features to the details of the local molecular structure. 
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